package scu.maqiang.heat;

import scu.maqiang.fes.BVPType;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.MVO;

import java.util.Arrays;

public class EarthData {
    public static void main(String[] args) {
        //heatTransfer();
        perFunc();
    }

    public static void perFunc() {
        String fileName = "O3w_depth_hor.txt";
        double[][] data = MVO.readArray2(fileName);
        int n = data.length;
        double[][] nodes = new double[n][];
        double[] height = new double[n];
        for(int i = 0; i < n; i++) {
            nodes[i] = new double[]{data[i][0], data[i][1]};
            height[i] = data[i][2];
        }
        System.out.println(Arrays.toString(data[11078]));
        int xNum = (int)Math.sqrt(n) - 1;
        Mesh2T3 mesh = new Mesh2T3().fromRectangleNodes(nodes, xNum, xNum, T3Type.Left);
        mesh.displayElements();
        mesh.toTecplot("earthData.dat", height);

        Mesh2Q4 mesh2 = new Mesh2Q4().fromRectangleNodes(nodes, xNum, xNum);
        mesh2.toTecplot("earthData2.dat", height);

        ScalarFunc perbutation = (xy, label, param) -> {
            double x = nodes[11078][0];
            double y = nodes[11078][1];
            double r = Math.sqrt((xy[0] - x) * (xy[0] - x) + (xy[1] - y) * (xy[1] - y));
            return 100 * Math.exp(- r * r / 4500);
        };

        double[] perb = mesh2.functionInMesh(perbutation, null);
        Arrays.setAll(height, i -> height[i] + perb[i]);
        mesh.toTecplot("perbutationData.dat", height);
    }
    public static void heatTransfer() {
        String fileName = "O3w_depth_hor.txt";
        double[][] data = MVO.readArray2(fileName);
        int n = data.length;
        double[][] nodes = new double[n][];
        double[] height = new double[n];
        for(int i = 0; i < n; i++) {
            nodes[i] = new double[]{data[i][0], data[i][1]};
            height[i] = data[i][2];
        }
        int xNum = (int)Math.sqrt(n) - 1;
        Mesh2T3 mesh = new Mesh2T3().fromRectangleNodes(nodes, xNum, xNum, T3Type.Left);
        mesh.toTecplot("earthData.dat", height);
        System.out.println(Arrays.toString(nodes[11078]));

//        Mesh2Q4Generator mg2 = new Mesh2Q4Generator();
//        Mesh2DQ4 mesh2 = mg2.fromRectangleNodes(nodes, xNum, xNum);
//        mesh2.toTecplot("earthData2.dat", height);


//        FES2T31 fs = new FES2T31(mesh);
//        SRMatrix A = new SRMatrix(fs.GetNdof());
//        SRMatrix M = new SRMatrix(fs.GetNdof());
//        fs.assembleHeatStiff(new double[]{1.0}, BVPType.COMMON, A);
//        fs.assembleHeatMass(new double[]{1.0}, BVPType.CONSIST_MATRIX, A);
        mesh.scale(1.0/4000, 1.0/4000);

        Heat2DT31 heat = new Heat2DT31(mesh);
        ScalarFunc source = (xy, label, param) -> 1.0 ;
        ScalarFunc zeroFunc = (xy, label, param) -> 0.0;
        ScalarFunc constOneFunc = (xy, label, param) -> 1.0;
        double[] x0 = mesh.functionInMesh(zeroFunc, null);

        double theta = 0.5;
        double T = 2.0;
        double dt = T / 400.0;
        double[] paramSource = {0.0};
        double[] paramfixed = {0.0};
        heat.GeneralHeatProblem(new BVPType[] { BVPType.COMMON, BVPType.CONSIST_MASS},
                constOneFunc, null,
                constOneFunc, null,
                source, paramSource,
                new double[] { theta, T, dt },
                new int[] { 1, 2, 3, 4 }, zeroFunc, paramfixed,
                null, null, null,
                null, null, null,
                null, null, x0,
                "MyHeat2DT31.dat");
//        Arrays.setAll(height, i -> height[i] + x0[i]);
//        mesh.toTecplot("perbutationData3.dat", height);
    }

}
